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Abstract 

Large  PEM  cells  will  be  used  in  future  proton  exchange  membrane  fuel  cell  (PEMFC)  power  plants  and  appropriate  tools  are  therefore  be  needed 
to  study  their  behaviour.  One  approach  to  understanding  single  cell  behaviour  involves  using  mathematical  models.  The  numerous  techniques 
used  in  this  work  to  describe  PEM  electrode  behaviour  require  different  scientific  disciplines:  chemical  engineering  and  electrochemistry.  This 
study  proposes  combining  residence  time  distribution  (RTD)  and  electrochemical  impedance  spectroscopy  (EIS).  The  investigation  focuses  on 
cathodic  DC  and  AC  responses  where  over-voltage  is  critical.  Results  demonstrate  that  although  gas  distribution  does  not  cause  additional  loops 
on  impedance  diagrams,  it  is  strongly  related  to  both  the  shape  and  amplitude  of  these  diagrams.  The  simulations  have  drawn  attention  to  operating 
conditions  that  can  threaten  the  life  of  the  PEM  cell:  under  these  setting  points  EIS  method  is  not  sufficient  to  detect  this  risk. 
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1.  Introduction 

In  recent  years,  fuel  cell  technology  has  attracted  consid¬ 
erable  attention  from  several  fields  of  scientific  research.  Fuel 
cells  are  highly  efficient  in  terms  of  energy  produced,  they  emit 
little  noise,  and,  if  fed  by  hydrogen,  are  non-polluting.  Fuel 
cell  development  is  particularly  important  for  stationary  and 
portable  applications.  Amongst  the  several  types  of  fuel  cells, 
distinguished  by  their  electrolyte,  PEMFCs  are  often  proposed 
because  of  their  mild  operating  conditions  (50-80  °C  temper¬ 
ature,  1-3  atm  pressure).  PEMFCs  are  composed  of  a  polymer 
electrolyte  sandwiched  between  two  porous  backing  layers  to 
form  a  membrane  electrode  assembly  (MEA).  The  electrodes 
(the  active  layers)  are  inserted  between  the  electrolyte  and  the 
backing  layers.  One  of  the  difficulties  in  characterising  PEM 
electrodes  is  related  to  the  scaling-up  of  the  fuel  cell.  Thus,  a 
specific  approach  must  be  developed  for  large  PEM  electrodes.  A 
well-established  powerful  tool  for  investigating  the  mechanisms 
of  physico-chemical  and  electrochemical  reactions  occurring 
at  PEMFC  electrodes  is  the  electrochemical  impedance  spec¬ 
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troscopy  (EIS).  This  method  allows  one  to  characterize  an 
electrode  behaviour  under  various  conditions  of  flow  feeding 
and  applied  current  density.  The  purpose  of  this  study  is  to  com¬ 
bine  the  continuous  stirred  tank  reactors  (CSTR)  approach  with 
DC  and  AC  electrochemical  descriptions.  Indeed,  it  is  possible 
to  improve  IES  explanations. 

According  to  Bernardi  and  Verbruge  [1]  and  Jaouen  et  al. 
[2],  since  water  is  produced  at  the  cathode,  the  impact  of  water 
management  on  cathode  performance  is  significant  at  high  cur¬ 
rent  densities.  Nevertheless,  these  authors  do  not  model  mass 
transfer  in  the  whole  MEA.  Experimental  studies  by  Buchi 
and  Srinivasan  [3]  confirmed  the  importance  of  gas  hydration: 
ionic  transport  through  the  membrane  is  the  most  limiting  phe¬ 
nomenon  when  gas  humidification  is  low.  However,  Sena  et  al. 
[4]  showed  that  the  access  of  gases  to  catalytic  surfaces  becomes 
a  limiting  factor  when  the  membrane  is  sufficiently  hydrated. 
Costamagna  and  Srinivasan  [5,6]  and  Okada  et  al.  [7]  confirmed 
these  observations  using  steady  state  and  transient  modelling. 
Typically,  the  electrode  (active  layer)  is  sandwiched  between 
the  polymer  membrane  (Nafion®)  and  the  gas  diffusion  layer 
(GDL).  On  both  sides  of  the  active  layer,  the  mass  balance  in  the 
GDL  and  the  polymeric  membrane  is  controlled  by  the  gas  flux 
inside  the  gas  channel  and  the  water  flux  through  the  MEA.  Most 
models  focus  only  on  the  cathode  side  of  the  fuel  cell,  because  the 
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Nomenclature 

bi  Tafel  slope,  (V  dec-1) 

ch20  water  concentration,  (molm3) 

Cdi  double  layer  capacitance,  (F  m-2) 

dpore  pore  diameter,  (m) 

Dij  diffusion  coefficient  of  species  i  through  j, 

(m2s-1) 

D?ff  effective  diffusion  coefficient,  (m2  s-1) 

Dm  effective  diffusion  coefficient  of  water  in  the 
membrane,  (m2  s-1) 

Df  Knudsen  diffusion  coefficient,  (m2  s-1) 

EW  equivalent  weight,  (kg  mol- 1 ) 

F  Faraday’s  constant,  (C  mol- 1 ) 

i  local  current  density,  (A  m-2) 

if  Faradaic  current  density,  (Am-2) 

i0(,-)  exchange  current  density,  (Am-2) 

J  average  current  density,  (A  m-2) 

L  thickness,  (m) 

Mi  molar  weight,  (kg  mol- 1 ) 

N  density  of  molar  flux,  (mol  m-2  s) 

P  Pressure,  (Pa) 

Psat  saturated  vapour  pressure,  (Pa) 

ro2  stoichiometry  oxygen  coefficient  (-) 

R  universal  gas  constant,  (J  mol-1  K-1)  (-) 

S  electrode  geometric  area,  (m2) 

t  time  component,  (s) 

T  temperature,  (K) 

x  space  component,  (m) 

X  electrode  area  ratio  (-) 

yi  molar  rate  of  the  gas  species  i,  (-) 

Z  impedance  (Q  m2) 

Greek  symbols 
e  porosity,  (-) 

y  roughness  factor,  (-) 

Oi  overpotential,  (V) 

i -p  density  of  molar  flux,  (mol  m-2  s) 

X  water  content  in  the  membrane,  (-) 

v  stoichiometric  coefficient  (-) 

Pdry  dry  Nafion®  density,  (kgm-3) 

crH+  protonic  conductivity,  (S  m-1) 

rr®  +  effective  protonic  conductivity,  (S  m-1) 

x  tortuosity,  (-) 

x0  osmotic  coefficient,  (mol  mol- 1 ) 

f  liquid/gas  flux  ratio  (-) 

^Nation  agglomerate  polymer  content  (-) 

Indices  and  exposits 
a  Anode 

act  Active  layers 

b  Backing  layers 

c  Cathode 

H2  Flydrogen 

H20  Water 

k  Space  discrete  component 


m  Membrane 

n  Space  discrete  component 

N2  Nitrogen 

02  Oxygen 


cathode  activation  potential  is  the  single  largest  source  of  inef¬ 
ficiency  in  fuel  cells.  Oxygen-reduction  reaction  at  the  cathode 
(1)  is  a  key  to  PEMFC  electrochemical  behaviour: 

02  +  4H+  +  4e-  — >  H20  (1) 

The  literature  proposes  various  ways  of  predicting  the  elec¬ 
trochemical  behaviour  of  electrodes.  Thus,  Springer  et  al.  [8,9] 
proposed  analytical  formulations  of  EIS  and  emphasized  the 
capacitive  and  resistive  effects  of  oxygen  reduction,  and  Eik- 
erling  and  Komyshev  [10]  complemented  these  steady  state 
models  with  their  descriptions  of  electrode  microstructure.  Stud¬ 
ies  conducted  by  Bultel  et  al.  [11]  showed  that  migration  and 
diffusion  processes  play  a  crucial  role  in  the  active  layer.  Song 
et  al.  [12],  by  taking  into  account  the  multi-step  reactions,  have 
also  shown  the  diffusion  limitation  occurring  in  the  backing 
layer.  The  results  of  these  studies  [11,12]  indicate  that,  on  a 
complex  plane  plot  (impedance  diagram),  diffusion  limitation 
leads  to  a  semi-circle  at  low  frequencies,  while  ionic  ohmic 
drop  produces  a  straight  line  with  a  45°  angle  at  high  frequen¬ 
cies,  which  is  in  agreement  with  Lasia’s  predictions  [13,14], 
Costamagna  and  Srinivasan  [5,6]  and  Djilali  et  al.  [15]  have 
suggested  that  multidimensional  modelling  improves  electrode 
description,  since  gas  composition  and  temperature  vary  along 
the  gas  channels.  In  addition,  these  phenomena  become  critical 
for  large  cell  areas.  Under  these  conditions,  recent  investigations 
[16-19]  have  been  carried  out  using  computational  fluid  dynam¬ 
ics  software.  Flowever,  resolution  in  these  approaches  requires 
substantial  computation  time.  Another  simple  way  to  describe 
the  fluid  dynamics  is  to  use  residence  time  distribution  (RTD) 
techniques.  Boillot  et  al.  [20]  have  demonstrated  that  the  bipo¬ 
lar  plates  are  similar  to  a  series  of  CSTRs,  while  Benziger  et  al. 
[2 1  ]  have  used  the  stirred  tank  reactor  (STR)  approach  to  explore 
the  auto-humidification  of  the  PEMFC.  Furthermore,  the  RTD 
model  does  not  require  lengthy  computing  times. 

The  present  study  uses  series  of  CSTRs  based  on  RTD  results 
obtained  on  PEM  electrodes  [20]  and  therefore  the  RTD  tech¬ 
nique  is  not  investigated  here.  The  purpose  of  this  study  is  to 
combine  the  CSTR  approach  with  more  classical  DC  and  AC 
electrochemical  descriptions.  AC  behaviour  is  currently  used  to 
characterize  the  electrochemical  generator:  the  electrochemical 
impedance  spectroscopy  technique.  It  is  therefore  important  to 
accurately  evaluate  EIS  responses.  Thus,  the  aim  of  this  study  is 
to  develop  an  approach  to  describe  the  behaviour  of  electrodes 
with  a  large  surface  area.  This  investigation  takes  into  account 
Stefan-Maxwell  diffusion  process  through  the  backing  and  cat¬ 
alyst  layers  [8]  completed  by  a  supplementary  term  for  Knudsen 
diffusion  in  the  small  pores.  In  an  earlier  study,  Ramousse  et  al. 
[22]  presented  a  theoretical  one-dimensional  model  in  steady 
state.  To  complete  this  description  of  the  PEM  electrode,  as  well 
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as  that  of  the  processes  occurring  therein,  the  present  study  pro¬ 
poses  AC  simulations.  Moreover,  the  CTR  approach  presented 
here  is  applicable  to  one-  or  two-dimensional  simulations  [23]. 

2.  Theory 


and  the  global  mass  balance  for  each  CSTR  is  given  by: 

4’in  -  4'0Ut  =  Sk(<Po2  +  /h2o)  (5) 

where  NJ  =  y^N^,  and  the  gas  flux  of  consumption  and  produc¬ 
tion  are  expressed  as  follows: 


On  both  sides  of  the  MEA,  the  bipolar  plates  combine  three 
functions:  thermal  control  through  cooling  water  channels,  gas 
feeding  via  the  gas  channels,  and  current  collector.  The  porous 
backing  layers  ensure  a  homogeneous  gas  distribution  on  the 
MEA.  The  active  layers  are  thin  layers  of  catalyst  dispersed 
on  carbon  particles  embedded  in  the  polymer  electrolyte.  The 
catalysed  electrochemical  and  chemical  reactions  occur  in  these 
layers.  Only  protons  and  water  can  go  through  the  polymer 
membrane  located  at  the  centre  of  the  assembly. 

2.7.  Gas  channel  model 


^membrane,  fc 

1njh2o 


2  F 


(6) 


(7) 


The  present  work  is  conducted  on  the  gas  phase  only  and  we 
have  introduced  representing  the  ratio  of  gas  molar  flux  to 
liquid  molar  flux,  which  is  estimated  using  a  numerical  optimi¬ 
sation  (1  <  %k  >  0).  The  saturated  vapour  pressure  of  water  is 
expressed  as  follows: 
psat  _  ^(23. 1964— (3816.44/  T— 46.13)) 


The  mass  balance  inside  the  gas  channel  is  solved  using  a 
simplified  fluid  dynamic  representation  such  as  a  CSTR  [20] 
series.  The  electrode  is  divided  into  k  elementary  units.  With 
such  a  hypothesis,  gas  transport  in  the  gas  channel  is  assumed 
to  be  infinitely  fast,  leading  to  homogeneous  gas  composition 
for  each  CSTR.  In  particular,  the  partial  pressure  of  a  compo¬ 
nent  in  the  elementary  unit  ( k )  is  equal  to  that  of  the  outlet 
and  to  that  of  the  elementary  unit  (k+1)  inlet,  as  illustrated  in 

Fig.  1. 

All  elementary  units  have  the  same  voltage  U  due  to  the 
high  electrical  conductivity  of  bipolar  plates  and  backing  lay¬ 
ers  (GDL).  This  implies  a  current  density  distribution  along 
the  gas  channel  that  is  determined  by  local  characteristics  J  — 
f(U,  Po2,  1%,  /jh2o,  T).  Thus,  the  partial  pressures  of  each 
component  in  each  elementary  unit  vary  and  make  it  necessary 
to  solve  mass  balances  and  current  densities  in  each  elementary 
unit  at  the  same  time.  The  current  voltage  expression  used  is 
that  proposed  by  our  electrode  model  [22],  The  molar  balance 
equation  of  each  CSTR  with  an  electrode  geometric  area  of  Sk 
for  each  component  in  the  gas  phase  in  steady  state,  is  written 
as: 


nh2o  -  nh2o 

—  £asVh2o 

(2) 

N^-Ngut 

II 

-6*" 

(3) 

=  0 

(4) 

2.2.  Electrode  model 

Srinivansan  and  Hurwirtz  [24]  assert  that  it  is  preferable  to 
use  a  volumetric  description  of  the  active  layer  and  a  formalism 
of  the  Butler- Volmer  type  rather  than  a  description  derived  from 
Tafel’s  law.  The  agglomerate  model  [25]  that  assumes  the  pres¬ 
ence  of  macro-pores  in  the  electrodes  seems  to  be  well  adapted 
to  the  PEMFC  anode  and  cathode.  Siegel  et  al.  [26]  used  it  to 
simulate  the  electrochemical  behaviour  of  a  MEA. 

Fig.  2  shows  a  schematic  representation  of  the  simplified  pore 
structure  assumed  in  the  agglomerate  model.  The  solid  phase 
consists  of  a  homogeneous  mixture  of  polymer  electrolyte  (e.g. 
Nafion®),  carbon,  and  catalyst  (Pt).  In  the  agglomerate,  there 
is  both  an  ionic  current  density  flux  and  an  electronic  current 
density  flux.  These  vary  with  the  x  coordinate.  The  ionic  current 
density  is  nil  at  the  active  layer/backing  layer  interface  (x  =  0) 
and  the  electronic  density  is  nil  at  the  membrane/active  layer 
interface  (x  =  Lact).  In  steady  state,  the  local  variation  in  elec¬ 
tronic  current  density  is  the  mathematical  opposite  of  the  local 
variation  in  ionic  current  density;  it  is  also  equal  to  the  faradaic 
current  density  if.  Cell  current  density  J  is  the  sum  (over  the 
electrode  thickness  Lact)  of  faradic  current  density: 


According  to  Kim  et  al.  [27]  the  oxygen-reduction  process 
is  limiting.  Boyer  et  al.  [28]  claim  that  the  ionic  conductivity 
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Fig.  2.  Schematic  description  of  the  PEMFC  half  cell  (cathode). 


of  the  active  layer  influences  the  electrochemical  performance 
of  the  electrodes.  However,  oxygen  access  to  the  catalyst  sites 
remains  problematic,  not  only  because  of  the  flooding  risk, 
but  also  due  to  its  low  partial  pressure.  As  Bernardi  and  Ver- 
bruge  [1]  proposed,  a  one-dimensional  model  can  describe 
cathode  operation  by  taking  into  account  only  the  oxygen 
diffusion  in  the  gas  pores  and  the  reduction  kinetics.  Under 
this  condition,  our  main  assumptions  in  the  model  are  the 
following: 

1 .  Mass  transfer  of  H2O  and  CL  is  considered  only  in  the  gas 
pores;  there  is  no  diffusion  of  water  in  the  agglomerate. 

2.  H2O  concentration  in  the  agglomerate  is  imposed  by  ther¬ 
modynamic  equilibrium  with  the  gas  phase. 

3.  The  ionic  conductivity  of  the  agglomerate  is  a  function  of 
the  ionic  conductivity  of  the  polymer,  which  depends  on  its 
water  content. 

4.  Only  the  migration  of  protons  is  considered  in  the  agglom¬ 
erate. 

5.  Pressure  and  temperature  are  uniform. 


[22]: 

<7pore  =  S  (12) 

1  -  £act 


Given  that  the  inside  of  the  backing  layers  is  porous,  it  is 
necessary  to  use  effective  diffusion  coefficients  [22]: 

Dfj  =  sl/3Dij  (13) 

Dq,  the  coefficient  diffusions  in  the  gas  phase,  come  from 
[29] 

Let  us  add:  YleJo  =  '  and  the  condition  that  only  water  can 
go  through  the  membrane,  in  a  direction  depending  on  current 
density  and  gas  hydration  [7],  These  are  given  by  the  following 
expressions: 


3  (pg  P  dye 

~dx  ~  RT~di 


Oxygen  consumption  and  water  production  in  the  porous 
active  layer  yield  the  following  mass  balance  equations: 

=  if_ 

dx  RT  dt  nF 


where  v  equals  to  0  for  nitrogen,  —  1  for  oxygen  and  + 1  for  water. 

The  overpotential  distribution  Eq.  (18)  in  the  active  layer  can 
be  calculated  by  means  of  the  proton  balance  (16)  and  Ohm’s 
law  (17)  in  the  ionic  phase: 


3(4+317*)  y  ,k  y  „  3 r,k 

3x2  -  f  -  z^Cdlir 


(16) 

(17) 

(18) 


The  effective  ionic  conductivity  of  the  agglomerate  depends 
on  Nafion® ’s  conductivity  and  is  simply  corrected  as  a  function 
of  the  agglomerate  polymer  content,  as  follows: 

4+M  =  VfNafionO-H+(x)  (19) 


Stefan-Maxwell  equations  [8],  completed  by  a  supplemen¬ 
tary  term  for  Knudsen  diffusion  in  the  small  pores  of  the  active 
layer  and  backing  layer  [22],  describe  the  diffusion  of  hydrogen 
and  water  in  the  gas  phase: 


dy*  _  RT  ( v  [  yke4  -  M 1  4  \ 

**  p  Ir  r  ^  1  D$) 


(10) 


The  Knudsen  diffusion  coefficient  Df  is  expressed  as  fol¬ 
lows: 


Tlf  —  c/noi 


(11) 


The  ionic  conductivity  of  the  membrane  depends  on  its  tem¬ 
perature  and  water  content,  according  to  this  correlation  [22]: 

<rH+  =  e(_£A(1/r“1/353))(0.00  13A.3  +  0.02987. 2  +  0.26587) 

(20) 

with 

Ea  =  2640e(“0<y-)  +1183  (21) 

where  the  polymer  water  content  7  is  defined  as  the  number  of 
water  moles  per  mole  of  sulfonic  acid: 

EW 

7  =  - cH2o  (22) 

Pdry 


where  eact  is  the  porosity  and  r  the  tortuosity  of  the  active  where  EW  (equivalent  weight)  represents  the  dry  membrane 
layer.  The  pore  equivalent  diameter  (<fpore)  can  be  estimated  by  weight  per  mole  of  the  sulfonate  group  (kg  mol-1);  p<jry  is  the 
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density  of  dry  polymer  (kgm-3).  The  water  content  (A.)  of  the 
polymer  electrolyte/gas  interface  assumes  thermodynamic  equi¬ 
librium  between  water  vapour  in  the  electrodes  and  liquid  water 
in  the  polymer.  The  sorption  curve  of  Hinatsu  et  al.  [29]  is: 

A.  =  0.3  +  10.8  -  16(^°)2  +  14.1 

(24) 


Thus,  the  agglomerate  ionic  conductivity  is  a  function  of  the 
x  co-ordinate.  Finally,  the  Butler- Volmer  law  [22]  is  generally 
used  to  describe  oxygen-reduction  kinetics  in  acid  media  (e.g. 
Nafion®): 


^(^Lback+Lact)  -  °> 

molar  flux  density  of  water  is  given  by  the  solution  of  the  Eq. 
k  membrane,  k 

(28)  ^H20(r=/.back+Lact)  -  <Ph20 


In  steady  state,  the  sum  of  diffusion  and  electro-osmotic 
fluxes  (Eq.  (26))  yields  the  following  first-order  differential 
equation: 


membrane,  k 

^h2o 


Pdry  dA* 

EW  dxT 


The  solution  of  Eq.  (27)  is: 


membrane,  k 

^h2o 


F 


r  4-a*  i 

[  a  1  —  exp(k^Lm)J 


(27) 


(28) 


4  =  io(02)  -  e(" 2.3/*Vl)j  (25) 

where  y0z  is  the  hydrogen  concentration  at  the  inlet  under 
open-circuit  conditions. 


2.3.  Membrane  model 


The  phenomenological  model  of  water  transport  in  the  mem¬ 
brane  proposed  by  Springer  et  al.  [8]  is  used  here.  According  to 
Okada  et  al.  [7] ,  the  water  electro-osmotic  flux  through  the  mem¬ 
brane,  which  is  always  directed  from  anode  to  cathode,  is  a  linear 
function  of  the  proton  flux  imposed  by  current  density:  The  water 
diffusion  flux  through  the  membrane  depends  linearly  on  the 
water  concentration  gradient.  The  sum  of  diffusion  and  electro- 
osmotic  fluxes  yields  a  differential  equation  for  the  A  variable: 


a2A* 


ik  dkk 

o— -z — i 

F  dx 


dXk 

Ik 


(26) 


where  to  is  the  electro-osmotic  drag  coefficient  and  Dm  the 
back  diffusion  coefficient  [7], 


2.4.  DC  solutions 

The  mass  balance  of  CSTR  series  uses  the  following  initial 
conditions:  yo2,N2  h20>  and  the  results  of  the  microscopic  model 
(ik,  N"anea).  The  steady  state  equations  of  each  elementary 
unit  are  solved  using  the  simplex  numerical  optimisation 
method. 

The  resolutions  of  mass  balances  are  numerically  obtained 
thanks  to  the  continuity  of  the  flux  between  the  backing  layer, 
the  active  layer  and  the  membrane.  The  boundary  conditions  of 
Eqs.  (14)  and  (15)  are  given  by  the  output  of  the  mass  balance 
in  CSTR: 


with 

hk  _  EWr0ik  ' 

Pdry  Dm  F 

Xk  and  Xk  stand  for  the  polymer  water  content  at  the  mem¬ 
brane/active  layer  interface.  Then,  assuming  thermodynamic 
equilibrium  between  water  vapour  in  the  backing  layers  and 
liquid  water  in  the  polymer,  the  sorption  curve  of  Hinatsu  et  al. 
[7],  Eq.  (24)  is  used  to  relate  partial  pressure  values  to  Xk  and 
A*. 

In  steady  state,  the  solution  of  Eq.  (18)  describes  the  distribu¬ 
tions  of  concentrations  and  electrode  overpotentials  as  functions 
of  x  coordinate.  This  differential  equation  has  to  be  solved  with 
the  following  boundary  conditions: 

-  Ionic  current  density  is  nil  at  the  active  layer/backing  layer 
interface  -rr^+d?)/dx|A=/Jtack  =  0. 

-  Electrode  overpotential  at  the  membrane/active  layer  interface 
is  known:  0(x=Lh!lcy+LM)  =  00  ■  In  practice,  r/0  is  adjusted  until 
J  equals  the  value  used  for  the  determination  of  hydrogen 
consumption  (29). 


This  set  of  equations  is  solved  by  iterative  methods  and 
only  needs  the  following  inlet  conditions:  y02,  yNz,  yH20 
and  Nn2°,N'^o,No20  where  N^)°  =  y^N™’0.  The  conventional 
way  of  investigating  the  influence  of  mass  flux  on  electrochem¬ 
ical  behaviour  is  by  using  the  inlet  molar  fluxes,  which  are 
proportional  to  the  current  density  and  the  stoichiometry,  for 
the  cathode  such  that: 


tXQ2 

y°o2 


J 

4 F 


(29) 


where  J  is  the  average  current  density  of  ik  and  ST  the  total 
geometric  area  of  electrode  ST  =  In  this  study,  the  initial 

inlet  molar  flux  density  Nq’°  is  fitted  to  the  case  where  there  is 
only  one  CSTR. 


concentration  of  each  species  at  the  gas  channel/backing  layer 
interface, 

molar  flux  density  of  oxygen  and  nitrogen  are  nil 
at  membrane/electrode  layer  interface  <PN2(x=z,back+Lact)  = 


2.5.  AC  solutions 

When  a  small  sinusoidal  amplitude  modulation  A  i  is  added 
to  the  applied  current  density  i,  the  concentrations  of  species 
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in  solution  or  adsorbed  are  also  modulated  in  a  sinusoidal 
way  around  their  steady  state  values.  Under  this  condition,  the 
faradaic  current  density  could  be  expressed  by: 


The  variation  of  oxygen  consumption  and  water  production 
in  the  porous  active  layer  yield  the  following  mass  balance 
equations: 


Aikf  =  ikf(t)  -  ik  =  d/fArfd )  +  dyo2  ikfAykQ2 (t)  (30) 


where  d^if,  dgif  et  3 yif  are  the  partial  derivatives  of  the  faradaic 
current  density  calculated  at  the  steady  state  such  that: 


V/  = 


dif 

dVc 


1 yo2 


ioah) 

bo2 


-2.3 


,2.ybo2\r,ci\yCh  +  hhe-2.3/b*\r,c 

y°o2  b* 


3  i  /&A|  =^2.37^1,01 

"°2  f  3y02  |,a  - 


(31) 


(32) 


We  denote  Laplace  transforms  (introducing  s  =  ja>  where  w 
is  the  angular  frequency)  of  the  time-dependent  perturbed  vari¬ 
ables  with  an  overbar. 


dx 


(39) 


These  ordinary  differential  equations  have  to  be  solved  in  the 
Laplace  domain  with  the  following  boundary  conditions: 


•  at  the  backing  layer/gas  channel  interface,  the  variation  in 
concentrations  is  nil  due  to  the  assumption  of  infinitely  fast 
mixing,  leading  to  an  homogeneous  gas  composition  (CSTR 
assumptions)  AykOi(x=0)  =  Ay£,2(?=0)  =  4y^O(x=0)  =  0, 

•  at  the  active  layer/membrane  interface,  the  flux  varia¬ 
tion  is  nil  for  the  oxygen  and  n i trogen  A(j^2 = 

AV>o  2(*=Lback+Lact)  =  °> 

•  according  to  steady  state  boundary  conditions,  the  flux  vari¬ 
ation  of  the  water  at  the  active  layer/membrane  interface  is 
determined  by  AVv\2o(x=wiCk+Lia)  = 


Atf  =  dyo2ikfAyk02  +  d/fAr,k  (33) 

Due  to  this  same  variation  in  the  cathodic  overpotential,  after 
Laplace  transforms  Ohm’s  law  can  be  written  for  any  point  of 
the  active  layer  as  follows: 


d Afjk  _  Aik 


(34) 


According  to  Okada  et  al.  [7],  the  variation  of  membrane 
water  content  in  the  Laplace  domain  involves  the  following 
expression  for  water  transport: 


d2AXk  ( ik  dAXk  3 Xk  Alk\ 

^  =  +  --j  (40) 

and  the  variation  of  water  transport  at  this  interface  results  in  the 
sorption  equilibrium  such  that: 


and  the  proton  balance  can  be  expressed  by: 
d  Aik  v  -i  i 

—  =  -Z-(A  ik  +  s  CW)  (35) 

To  solve  this  system  numerically  with  simple  discretization, 
we  use  a  constant  step  size  h.  The  discreet  solution  may  be 
written  as: 


—k - £ - =  K(T)  (41) 

Ayn2o(x=LbaA+Laa) 

In  order  to  determine  the  impedance  of  each  elementary  elec¬ 
trode  (Zk),  the  active  layer  was  divided  into  n  elements  (length 
h),  where  n  varies  from  500  to  1000.  Thus,  the  total  impedance 
is  given  by: 


l/Zk_1  +  yh(l/Z^+  scd i) 


azi ye  ki  /  ^ 

dx  =  ~P  ^ 


{f>\&%+y\Av\}-^Ay\  +  y#Ayke)] 


Z  SkA-JZk 


(42) 


where  Zkn  =  Arjk/Aik  is  the  local  impedance  at  n  active  layer 
abscise  and  =  Atjk/Ail ^  the  faradic  impedance  at  n.  The 
total  impedance  of  the  electrode  can  be  represented  in  the  form 
of  a  line  transfer  [19,10,11],  Under  these  same  conditions,  the 
mass  Eqs.  (8)-(14)  become: 


To  simulate  the  EIS ,  various  inlet  fluxes  are  investigated  under 
many  different  levels  of  relative  humidity.  But  the  most  impor¬ 
tant  contribution  of  this  study  is  the  description  of  fluid  dynamics 
made  possible  through  the  use  of  CSTR  models.  Different  archi¬ 
tectures  are  investigated  with  a  view  to  predicting  AC  electrode 
behaviour:  series  and  parallel  pathways  of  CSTR. 

3.  Results  and  discussion 


AVe\ 

»Pe) 


(37) 


This  equation  of  diffusion  transport  in  the  Laplace  domain 
confirms  that  J2gAyjg  =  0.  The  mass  balance  inside  the  GDL  is 
expressed  as  follows: 


d^<PN2 

dx 


ta^2 


(38) 


The  values  of  selected  parameters  used  to  simulate  both  DC 
and  AC  behaviour  are  listed  in  Table  1 .  Most  of  them  are  based  on 
data  from  the  literature.  The  range  selected  for  kinetic  constants 
was  chosen  in  order  to  obtain  cathodic  current  densities  simi¬ 
lar  to  those  observed  in  the  PEMFC  cathode  (for  \r]c\  =0.55  V, 
J<1  Am-2). 

In  this  study,  three  kinds  of  humidification  modes  are  inves¬ 
tigated:  auto-humidification,  humidification  of  inlet  gas  and 
humidification  via  the  membrane.  Moreover,  the  inlet  gas  flux 
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Table  1 

Simulation  values 


Parameter 

Value 

Refs. 

Lh»r* 

230  x  10“6  m 

[22] 

Lact 

10  x  10_6m 

[22] 

Tm 

125  x  10"6  m 

[22] 

£hnrk 

0.8 

[22] 

£act 

0.3 

[22] 

^Nation 

0.5 

[22] 

8 

3  x lO"6  m 

[26] 

S 

1.7- 

[22] 

y 

100  m2m"2 

[30] 

*02  =  *02 

120  mV  dec"1 

[30] 

10(02) 

4  x  10-7  Acm-2 

[30] 

Em 

0.28 

[22] 

EW 

1.1  kg  mol-1 

[7] 

P  dry 

2020  kg  m“3 

[7] 

TO 

2.5/22- 

[7] 

Dm 

3  x  10_10m2  s_1 

[7] 

Do2/H20 

3.20  x  10“5  m2  s-1 

[22] 

DoJ N2 

2.41  x  10“5  m2  s"1 

[22] 

Fig.  3.  Simulated  polarisation  resistance  as  a  function  of  ro2  for  dry  conditions, 
(a)  1  CSTR,  (b)  20  CSTR  with  (Q),  r0,  =  2.2;  (□),  r0%  =  \\. 


also  affects  the  electrochemical  behaviour  of  the  PEMFC  cath¬ 
ode.  This  flux  is  currently  expressed  using  the  stoichiometry 
oxygen  coefficient,  expressed  as  follows: 


r02  =  4  F 


(43) 


The  model  presented  here  does  not  take  into  account  the 
appearance  of  liquid  water.  Based  on  this  assumption,  the  values 
of  ro2  are  estimated  in  order  to  obtain  a  gaseous  condition  in  the 
set  of  calculated  cases.  A  simple  mass  balance  [31]  could  be 
used  to  estimate  the  minimum  flux  of  air  or  the  ro2  value  for  the 
non-appearance  of  liquid  water  on  the  cathode: 


r02(T)  =  0.21 


2+ P/P-  Pm(T) 

P/P  -  /Jsat(T)  -  RHqP/P  -  RHc/\at(T)J 

(44) 


The  limit  of  liquid  water  appearance  depends  on  the  relative 
humidity  (RHc)  of  inlet  gas.  This  mass  balance  equation  con¬ 
siders  that  the  water  produced  by  oxygen  reduction  is  totally 
evacuated  by  the  gas  flux.  This  relationship  does  not  take  into 
account  the  water  flux  provided  from  the  anode  side:  for  this 
reason,  a  high  stoichiometry  oxygen  coefficient  is  proposed. 
The  selected  values  of  ro2  =  2.2  and  1 1  should  avoid  flooding  at 
60  °C  for  any  simulated  cases  such  as  RHc  <0.8  and  RHc  <  0.98. 

3.1.  Electrochemical  behaviour  using  dry  inlet  conditions 

The  following  simulations  show  the  effects  of  two  phenom¬ 
ena:  first,  the  inlet  gas  flux  impact  and  second,  the  influence  of 
gas  distribution  on  the  electrode.  The  effects  of  gas  distribution 
are  investigated  using  the  CSTR  description:  for  each  CSTR, 
there  is  a  corresponding  uniform  electrode  zone.  In  Figs.  3  and  4, 


(a)  „  Irrl  =  0.25  V 

xio  * 


0  0.5  1  1.5  2 

real(Z)  /  Qm2  x  1°'3 


(b)  xio-4  1^1  =  0.25  V 
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. \ . 

>/ 

_ _ 1  : 
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real(Z)  /  Qm2  x  10"3 
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. 
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Fig.  4.  Simulated  complex  impedance  diagrams  as  a  function  of  ro2  for  dry  conditions  at  |ijc|  =  0.25V  and  |r/c|  =  0.55  V,  (a)  1  CSTR,  (b)  20  CSTR  with  (O), 
ro2  =  2.2;  (□),  ro2  =  11.  The  arrows  indicate  decreasing  frequency  values  from  100  kHz  to  10-3  Hz. 
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Fig.  5.  Simulation  of  local  impedance  diagrams  ro2  for  dry  conditions  with 
|  r/c  |  =  0.4  V  and  ro2  =  1 1  ■  The  arrow  indicates  decreasing  frequency  values  from 
100  kHz  to  10-3  Hz. 


the  simulation  of  polarisation  resistances  or  impedance  diagrams 
demonstrates  the  competition  between  proton  and  oxygen  access 
to  the  catalyst.  The  impedance  diagram  exhibits  the  shape  of 
finite  length  diffusion  impedance  when  proton  access  controls 
electrochemical  performance  [13]  and  that  of  a  semi-circle  when 
oxygen  access  controls  performance  [11].  In  the  case  of  a  cath¬ 
ode  with  operating  uniformity  (CSTR  =1),  the  limiting  transport 
step  is  oxygen  diffusion:  the  impedance  diagram  exhibits  a  semi¬ 
circle  shape  at  low  frequencies,  decreasing  with  the  increasing 
stoichiometry  oxygen  coefficient  (|/jc|  =0.25  V,  Fig.  4a).  The 
proton  limitation  could  appear  only  at  high  current  densities  with 
a  high  inlet  air  flux:  a  45°  straight  line  in  Fig.  4a  ( |  rjc  \  =  0.55  V). 
For  electrodes  with  operating  non-uniformities  (e.g.  20  CSTR), 
the  analysis  is  similar  to  low  cathodic  over-voltage  (|  rjc  \  =  0.25  V, 
Fig.  4b).  But  at  high  current  densities  (\rjc  \  =0.4  V,  Fig.  3b), 
proton  access  becomes  critical,  which  is  confirmed  by  increas¬ 
ing  electrode  resistance  as  the  stoichiometry  oxygen  coefficient 
increases  (Fig.  4b,  |/?c|  =0.55  V). 

A  study  of  local  impedance  is  necessary  to  understand  this 
behaviour.  Local  impedance  diagrams  are  shown  in  Figs.  5  and  6. 
The  CSTR  are  designated  as  1-20  from  the  inlet  gas  to  outlet 


Fig.  6.  Simulation  of  local  impedance  diagrams  ro2  for  dry  conditions  with 
\r]c\  =0.4  V  and  ro2  =  2.2,  The  arrow  indicates  decreasing  frequency  values 
from  100  kHz  to  10"3  Hz. 


gas.  These  figures  are  selected  to  highlight  ionic  transport  (a 
45°  straight  line)  and  oxygen  transport  (capacitive  loop  at  low 
frequencies).  In  both  figures,  the  impedance  diagrams  for  the  first 
CSTRs  reflect  the  dry  electrode  effects.  Thus,  these  zones  do  not 
produce  enough  water  to  satisfy  the  humidification  requirements 
of  the  catalytic  layer.  However,  subsequent  zones  receive  more 
humidity. 

Fig.  5  demonstrates  that  oxygen  access  does  not  constitute 
a  limitation,  and  that  electrochemical  performance  is  limited 
more  by  proton  access.  Under  these  conditions,  proton  access  is 
improved  along  the  gas  channel  (ro2  =  11)  and  local  impedance 
decreases  from  the  inlet  to  the  outlet,  gradually  becoming  con¬ 
stant.  These  simulation  results  are  similar  to  the  experimental 
observations  of  Schneider  et  al.  [32] ,  obtained  using  a  local  elec¬ 
trochemical  impedance  technique.  But  for  ro2  =  2.2  (Fig.  6), 
only  the  inlet  zone  exhibited  limitation  by  proton  access,  due 
to  the  dry  inlet  conditions.  The  other  zones  exhibited  limitation 
by  oxygen  access  (low-frequency  semi-circle).  In  addition,  this 
capacitive  contribution  increases  along  the  gas  channel  due  to 
the  consumption  of  oxygen. 

At  this  point  in  the  investigation,  it  is  important  to  understand 
the  influence  of  the  CSTR  description.  This  method  results  from 
the  dynamical  analysis  of  the  gas  distribution  in  the  bipolar  plates 
[20].  The  following  simulations  underline  the  dependence  of  AC 
responses  on  CSTR  modelling. 

The  simulations  shown  in  Fig.  7  demonstrate  that  the  impact 
of  flux  distribution  on  the  electrode  is  a  function  of  CSTR 
description.  In  this  figure,  two  cases  are  proposed:  low  inlet 
flux  of  air  (ro2  =  2.2)  and  high  inlet  flux  of  air  (ro2  =  11). 


Fig.  7.  Simulated  complex  impedance  diagrams  as  a  function  of  CSTR  descrip¬ 
tion  for  dry  conditions  with  (a)  ro2  =  2.2  and  (b)  ro2  =  11,  Itfcl  =0.45  V  (the 
arrow  indicates  an  increasing  number  of  CSTR  descriptions).  ( — ■),  1  CSTR; 
(-  -),  2  CSTR;  (-  -),  5  CSTR,  (. . .),  10  CSTR  and  (•),  20  CSTR.  The  arrow 
indicates  decreasing  frequency  values  from  100  kHz  to  10-3  Hz. 
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With  a  high  air  flux  (Fig.  7(b)),  the  inlet  zones  do  not  pro¬ 
duce  enough  water  to  ensure  a  correct  level  of  humidification: 
the  water  produced  is  moved  out  of  the  cell  too  quickly,  as 
shown  in  Fig.  5.  This  phenomenon  is  amplified  by  increasing 
the  number  of  homogenous  zones:  increasing  of  drying  zones. 
Thus,  proton  migration  into  the  active  layer  could  control  the 
oxygen-reduction  process,  as  observed  by  Jaouen  et  al.  [2,33]. 
On  the  other  hand,  at  low  inlet  flux  levels  the  limiting  process 
becomes  gas  access,  as  shown  in  Fig.  6.  Thus,  in  Fig.  7a,  the 
impedance  diagrams  decrease  as  the  number  of  homogenous 
zones  increases  due  to  better  oxygen  access  in  the  inlet  zones. 
In  these  conditions,  a  progressive  distribution  of  reactant  on 
electrodes  (20  CSTR)  improves  electrochemical  performance. 

3.2.  Influence  of  humidification  mode  on  EIS 

The  auto-humidification  of  the  electrode  was  studied  in  the 
Section  1  of  this  paper,  but  the  humidification  modes  of  PEM 
electrodes  also  have  an  impact  on  IES.  Three  types  of  humidifi¬ 
cation  are  investigated  in  this  study: 


frequencies  depends  on  the  inlet  gas  flux  and  on  humidification 
modes.  In  case  n°2  for  \rjc  \  =0.25  V,  the  value  of  the  apex  fre¬ 
quency  increases  from  14  Hz  (ro2  =  2.2)  to  22  Hz  (ro2  =  11), 
while  in  case  n°l  for  \i]c  \  =0.25  V,  this  value  increases  from 
17  Hz  (ro2  =  2.2)  to  29  Hz  (ro2  =  1 1).  This  observation  under¬ 
lines  the  competition  between  gas  access  (low  frequencies)  and 
proton  access  (high  frequencies). 

Hydration  of  the  cathode  via  the  anode  side  (case  n°3)  pro¬ 
duces  a  special  effect.  Under  this  condition,  the  electro-osmotic 
flux  is  added  to  the  back  diffusion  flux.  Oxygen  diffusion  through 
the  cathode  slows  down  with  the  increasing  water  flux  (Eq.  (37)), 
according  to  the  following  ratio: 

<4,0^2  or  o2  (45) 

DHio,N2  or  O2 

This  phenomenon  corresponds  to  an  increasing  semi-circle 
in  the  impedance  diagram.  In  general,  humidification  improves 
electrochemical  performance  and  minimizes  the  effect  of  the 
oxygen  stoichiometric  coefficient  (ro2  ). 


1.  auto-humidification,  extensively  discussed  in  Section  3.1 
where  there  is  both  a  dry  inlet  gas  and  dry  conditions  at  the 
anode/membrane  interface, 

2.  humidification  of  inlet  gas  at  half  the  saturated  water  pressure 
and  dry  conditions  at  the  anode/membrane  interface, 

3.  non-humidification  of  inlet  gas  but  the  anode/membrane 
interface  is  saturated  by  the  water. 

The  45°  straight  line  appears  at  high  frequencies  for  high  feed 
fluxes  and  high  current  densities  (Fig.  8).  However,  the  change 
in  shape  of  the  impedance  diagrams  is  faster  for  10  CSTR  than 
for  1  CSTR.  This  behaviour  is  linked  to  the  existence  of  dry  local 
effects  that  were  presented  in  Figs.  5  and  6.  Of  course,  the  apex 
frequencies  increase  with  over-voltage,  but.  the  value  of  these 


3.3.  Influence  of  flux  distribution  on  EIS 

In  this  section,  a  parallel  pathway  for  the  gas  flux  is  inves¬ 
tigated  (Fig.  9).  A  new  operating  condition  is  proposed:  the 
temperature  is  set  at  65  °C  to  ensure  total  consumption  of  the 
oxygen  and  better  steam  transport  (  sat”  w  0-21  atm).  To 
understand  the  effects  of  gas  distribution  on  the  PEM  electrode 
(case  of  parallel  pathways),  an  additional  parameter  is  added 
(Eq.  (46)).  In  an  earlier  section  of  this  work  (Section  3.1),  all 
CSTRs  had  the  same  electrode  area,  while  in  this  section,  the 
area  of  each  CSTR  is  deduced  from  the  following  expression: 

4  =  yx  100  (46) 


Fig.  8.  Simulated  complex  impedance  diagrams  as  a  function  of  humidification  conditions  8  for  |  ijc  |  =  0.45  V,  with  ro2  —  2.2  and  1 1 ;  (-□-),  humidification  of  inlet 
gas  at  1/2  of  the  saturated  water  pressure;  (B-  -),  water-saturated  anode;  (0)>  auto-humidification  (dry  conditions).  The  arrow  indicates  decreasing  frequency  values 
from  100  kHz  to  10-3Hz. 
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Fig.  9.  Molar  flux  distribution  in  elementary  units  (CSTR). 


In  this  section,  we  consider  two  CSTRs  such  as  1  =  + 

X2S  to  simulate  a  threat  to  the  life  of  the  PEM  cell  due  to  gas 
distribution.  Thus,  the  simulation  conditions  are  as  follows: 


4.  Conclusion 


•  Tq2  =  1 1  (principal  gas  pathway), 

•  /q  1  such  as  yQ2  ~  10-5  (case  of  electrode  dead  zones), 

•  on  both  CSTRs,  the  inlet  gas  is  dry, 

•  the  anode  side  of  CSTR  n°  1  is  saturated  by  the  water  (case  of 
water  accumulation  in  dead  end  mode  of  fuel  feed). 

Fig.  10  shows  the  effect  of  current  density  distribution  on 
polarisation  resistances.  As  expected,  polarisation  resistance 
increases  with  increasing  X2S  area.  Although  the  most  important 
part  of  the  current  is  provided  by  CSTR  n°  1,  CSTR  n°2  affects 
not  only  the  value  of  polarisation  resistance  but  also  the  shape  of 
the  impedance  diagrams  (Fig.  1 1):  in  Fig.  1 1  the  45°  straight  line 
disappears  for  X]s  =  60%.  Under  these  conditions,  it  is  difficult 
to  observe  local  drying  using  the  EIS  technique.  On  the  other 
hand,  when  liquid  water  is  accumulated  inside  the  electrode,  the 
limiting  diffusion  could  result  in  an  additional  low  frequency 
contribution  being  observed  in  the  impedance  diagram:  flooded 
electrode  [14],  or  low  electrode  porosity  [23] .  Thus,  it  is  possible 
that  hazardous  operating  conditions  cannot  always  be  estimated. 


Simulations  are  conducted  only  in  the  gas  phase  with  a  view 
to  determining  the  impact  of  gas  distribution  in  the  gas  channel 
on  impedance  diagrams.  This  work  has  presented  the  theoreti¬ 
cal  bases  for  an  experimental  method  for  characterising  large 
PEM  electrodes  using  the  RTD  method  (CSTR  model)  cou¬ 
pled  with  IES  techniques  (AC  model).  Thus,  to  improve  IES 
explanations,  it  is  necessary  to  understand  gas  distribution.  This 
approach  does  not  require  lengthy  computations  and  the  simula¬ 
tions  are  easier  to  fit  to  experimental  data.  In  all  the  simulations, 
only  one  capacitive  contribution  is  observed.  Theoretically,  AC 
active  layer  behaviour  involves  only  one  capacitive  loop  and  an 
additional  one  due  gas  access  through  backing  layer  could  not 
appear  (large  porosity  e  =  0.8).  Fluid  distribution  does  not  pro¬ 
vide  an  additional  loop  but  the  changing  shape  is  strongly  linked 
to  the  CSTR  description.  So,  it  is  difficult  to  observe  local  drying 
or  flooding  using  the  EIS  technique  alone.  Moreover,  operating 
conditions  that  can  threaten  the  life  of  the  PEM  cell  are  not  easy 
to  detect  using  only  the  EIS  method.  Thus,  this  work  has  drawn 
attention  to  the  relative  advantages  of  using  a  fluid  dynamics 
description  coupled  with  the  EIS  technique.  In  the  future,  it  will 
be  important  to  include  the  flooding  or  partial  flooding  effects 
that  were  not  taken  into  account  in  this  study. 
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